The Debye- Waller factor of liquid silica: Theory and simulation 
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We show that the prediction of mode-coupling theory for a 
model of a network-forming strong glass-former correctly de- 
scribes the wave- vector dependence of the Debye-Waller fac- 
tor. To obtain a good description it is important to take into 
account the triplet correlation function C3, which we evaluate 
from a computer simulation. Our results support the possi- 
bility that this theory is able to accurately describe the non- 
ergodicity parameters of simple as well as of network-forming 
liquids. 

PACS numbers: 61.20.Lc, 61.20.Ja, 64.70.Pf, 51.10.+y 

The quantitative description of the glassy dynamics 
in liquids is one important goal of modern research in 
condensed matter. Work in the last decade has pro- 
vided evidence Q that the so-called mode-coupling the- 
ory (MCT) is able to describe the slow dynamics of 
fragile liquids in the weakly supercooled state. Detailed 
theoretical predictions for model systems — including 
hard sphere systems, simple binary liquids, and molecu- 
lar liquids — have been found to be in remarkable agree- 
ment with experimental measurements as well as with 
simulation results l^-Q. E.g. the full ^-dependence of 
the Lamb-Mossbauer and of the Debye-Waller factors is 
predicted well by the theory. 

Recently, it has been shown that also intermedi- 
ate and strong glass-formers, such as glycerol or silica 
(S1O2) show features that are in qualitative agree- 

ment with the predictions of MCT. However, in the case 
of silica, a detailed comparison between theory and nu- 
merical data has questioned the ability of MCT to de- 
scribe correctly the Debye-Waller factors of this impor- 
tant network-forming material 1 1 O]. In this Letter we 
show that the disagreement between MD-results and the- 
oretical predictions — for the case of silica — were not 
due to failure of MCT to describe caging in this network 
forming liquid but to a further approximation which is 
assumed for the sake of simplicity in the commonly-used 
MCT equations. We show that, once this approxima- 
tion is avoided, MCT is able to describe accurately the 
cages in silica, opening the way for a full description of 
dynamics in network forming liquids above the critical 
temperature of MCT. 

To start we briefly review the MCT equations. The 
central quantity is the coherent intermediate scatter- 
ing function F(q, t), which for binary system can be 
written as a 2 x 2 matrix with entries Fij (q, i) = 



(<5pi(q, 0)<5p*(q, t)). Here Pi(q, t) is the density fluctua- 
tions for wave- vector q, and i is the label for the species. 
The equation of motion for F(q, t) is given by 



F(q,t)+n'(g)F(q,t)+ / drM(q, t- r)F(q, r) = 0, (1) 
Jo 

where the frequency matrix is given by \H 2 (qj\ . . = 

q 2 kBT(xi/m i )Y,k$ik[S~ 1 {q)\ k y Here q = |q|, Xi is 
the concentration of species i, mi is their mass, and 
S is the partial structure factor matrix. Within the 
framework of MCT, the memory function M is given 
by Mij — XiksTNij /rrii, where the matrix Nij is a 
quadratic form in Fij(q, t): 

3 J K 1 a/3 a'P' 

Vja/ft, (q, q - k)F aa , (k, t)F pp , (q - k, t). (2) 

Here n is the particle density, and the vertices (q, k) 
are given by 



q • (q - k) 



S la cf{q - k) 



+n<?C3 Q/3 (q,q-k) ■ Xi . (3) 



The function ^(q) = 5ij/xi — [S 1 (q)] ij is the direct 

correlation function, and C3 (q, k) is the triple correla- 
tion function, which is related to the triple density fluc- 
tuations via 

(Sp a (q)% (k)S P ; (q + k)) = N Sae (q) S fja {k) 

ear] 

x5 7I) (|q + k\)(5 ea S ev 5 na /x 2 e + n 2 4 ar > (q, k)). (4) 

The only input required by MCT are the static quanti- 
ties S, C3, and n. Temperature enters the equations only 
through the explicit T-dependence of fl 2 (q) and the im- 
plicit T-dependence of the static quantities. The solution 
of this type of equation shows at low T a two-step relax- 
ation dynamics, if T is above T c , the so-called critical 
temperature of MCT O . Below T c the time correlation 
functions do not decay to zero anymore, thus signaling 
that the system is not longer ergodic. Hence the height 
of the plateau in a time correlation function is usually 
called the non-ergodicity parameter (NEP) since it mea- 
sures that fraction of the correlation that does not decay 
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to zero even at long times. The physical relevance of 
the NEP can be inferred from the fact that it can be 
directly measured in light and neutron scattering exper- 
iments. For the case of F(q, t) the NEP is F c (g), the 
Dcbye- Waller factor for wave- vector q. 

In one of the first attempts to solve the MCT equa- 
tions for soft spheres, the possibility of setting in Eq. (|^) 
C3 = was considered 0. This approximation corre- 
sponds to a factorization of the triple density correlation 
in g-space in product of the three pair density correlation. 
It was found that for this simple liquid, this approxima- 
tion does not significantly affect the MCT predictions. 
All further works have build on this information. In 
this Letter we check the validity of this approximation 
for network- forming liquids. We use molecular dynam- 
ics (MD) simulations to determine the triple correlation 
function for silica Jl2| and use these functions to solve 
the full MCT equations. This calculation allows us to 
determine for the first time whether or not MCT is able 
to give a quantitative description of the cages of strong 
glass-formers. For completeness, we perform the same 
study for a well studied binary mixture of Lennard-Jones 
particles (BMLJ) flf|. 

In the past, the dynamics of both systems has been 
carefully analyzed and it has been shown that the relax- 
ation dynamics shows the qualitative features predicted 
by the MCT " 
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For the case of the BMLJ also 
a quantitative comparison has been made in that F c (q) 
has been calculated theoretically (assuming C3 = 0) and 
compared with simulation data. This comparison showed 
that very good agreement between theory and simulation 
is obtained, even if C3 = is assumed, in agreement with 
the conclusions of Barrat and Latz [pj . But for the case 
of Si02 a similar comparison showed that the MCT pre- 
diction for F c (a) does not describe properly the confining 
cage if c 3 = pjn). 

In order to investigate the influence of C3 it is nec- 
essary to determine c^ 71 (q, k) with high precision. For 
this we calculated {Sp a (k)Sp/3(p)Sp*(q}) with q = k + p 
for all triplets of values |k|,|p|,|q|, with |p| > |k| and 
p — k < |q| < p + k|, and then used Eq. (jj) to de- 
termine C3. We have selected a mesh size A, where A 
was 0.154A- 1 and 0.334 for the case of Si0 2 and BMLJ, 
respectively. Only wavevectors with modulus less than 
100A have been considered. For each triplet of value 
for |k|,|p|, |q| we randomly selected up to 10 2 combina- 
tions of wave-vectors k, p, q, calculated for these the 
density fluctuations and then the triple density correla- 
tion function. For BMLJ this was done for 12309 com- 
pletely independent configurations, corresponding to 100 
million time steps in the simulation. For SiC>2 we ana- 
lyzed 8577 configurations (20 million time steps). Note 
that this very large number of independent configura- 
tions was needed in order to determine C3 with suffi- 
cient high accuracy to make sure that the final results 
for the NEP do not depend anymore on the noise in C3. 



We mention that using only 1000 independent configu- 
rations was, e.g., not sufficient to guarantee this. Due 
to the large effort needed to generate the configurations 
and to analyze them, we determined C3 only for one tem- 
perature, namely T = 1.0 (BMLJ) and T = 4000K Si0 2 . 
Thus in the calculation of the MCT vertices in Eq. (|^) 
we assumed that the direct correlation functions depend 
on temperature but that the T-dependence of C3 can be 
neglected. 
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FIG. 1. Wave- vector dependence for the non-ergodicity 
parameters for the BMLJ system. The solid and dashed 
curves are the theoretical prediction with and without the 
inclusion of the C3 terms in Eq. t^. The circles are the MD 
results from 111811 - 
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FIG. 2. The same quantities as in Fig. |y, but for silica. The 
MD data is from Ref. [uj|. 

From the equations of motion, Eqs. (|l|)-@, we deter- 
mined T c and F c (q), using the iteration procedure from 
Ref. J||. The results obtained are presented in Figs. |l| 
and |2 were we show the three components of the NEP 
for the case of the BMLJ and silica, respectively. In 
each panel three curves are shown: the data from the 
simulation (filled circles) jl8), the theoretical prediction 
if C3 = (curves labeled with "02"), and the theoretical 
prediction if C3 is taken into account (curves labeled with 
"03") [^9[ . From Fig. |]we see that the curves for C2 and 
C3 are very close together, from which we conclude that 
in the case of the BMLJ the inclusion of C3 changes the 
prediction of the theory only weakly. The main difference 
between the two curves is in the BB correlation, since the 
C3-curve describes the simulation data even better than 
the C2-curve. Also the theoretical prediction for T c is es- 



sentially independent of whether or not C3 is neglected, 
in that T c decreases from T c = 0.922 to T c = 0.910 if 
C3 is taken into account. These results are in agreement 
with the findings of Ref. || . 

From Fig. ||we recognize that for the case of SiC>2 the 
influence of C3 on the NEP is much stronger than for the 
BMLJ system. For this network-forming liquid the inclu- 
sion of c 3 into the vertices leads to predicted F c (q) which 
are in substantially better agreement with the ones from 
the simulation than the ones from the C2-theory. The 
improvement is particularly noticeable for q above the 
location of the first peak, a length scale that corresponds 
to the distance between two neighboring tetrahedra. E.g. 
in the case of the Si-Si correlation the relative difference 
between the theoretical curve and the MD data is de- 
creased by a factor around five. The remaining differ- 
ence is now only on the order of the error of the MD 
data. Note that F c (q) provides information on the cage 
in which a particles is confined and which is formed by 
its neighboring particles. The fact that in the vicinity of 
the first maximum the two theoretical curves as well as 
the MD data are very close together shows that already 
the C2-theory is able to capture the structure of the cage 
on this length scale. For larger wave- vectors the C2-curve 
is too high, which means that the size of the cage is un- 
derestimated. Only if the terms due to C3 are taken into 
account, a reliable description for the cage is obtained. 

For the Si-Si correlation the mentioned decrease of the 
NEP is more pronounced than for the two other correla- 
tion functions, which is reasonable since it is the Si-atoms 
that sit in the center of the tetrahedra, i.e. that make up 
the network structure, and hence it can be expected that 
the inclusion of the C3-terms is important. Nevertheless, 
from the figure we see that the inclusion of these terms 
lead also to a significant improvement of the theoretical 
prediction for the NEPs for Si-0 and 0-0 if q is larger 
than the location of the first peak. E.g. for the case 
of 0-0 the relative error of the eveurve for q > 2.0A is 
only half as large as the one for the C2-curve. Thus we 
conclude that the C2-MCT is able to give a good descrip- 
tion of the cage for the oxygens only on the length scale 
corresponding to the first maximum of F(q, 0). 

The strong influence of the C3-terms on F c (q) is also 
found in the value of T c . If one sets C3 = one finds T c — 
3962K, whereas the C3-theory predicts T c = 4676K. Thus, 
whereas for the BMLJ system we found only a change 
of the order of 1%, we now find a change around 18%. 
We mention, however, that MCT usually overestimates 
the value of T c , and indeed we have T^ 10 — 0.435 and 
T* ID = 3330K for BMLJ and Si0 2 , respectively Jl7],||. 
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FIG. 3. Wave-vector dependence of the A-A and Si-Si com- 
ponent of the memory function (upper and lower panel, re- 
spectively). The dashed and solid curves are for the version 
of MCT in which the C3-terms are neglected and taken into 
account, respectively. 

Finally we show in Fig. || the g-dependence of Naa 
and NgiSi from Eq. (||) for the corresponding T c . From 
this figure we see that for the BMLJ case the C2-memory 
function is very similar to the one in which C3 is taken 
into account and thus it is not surprising that also the 
NEP and the T c are not changing. This is in stark con- 
trast to the SiC>2 system, for which the memory func- 
tions from the c%- and C3-theory are very different. Note 
that in a binary mixture each component of the memory 
function has contributions from 16 different terms (see 
Eq. (||)). We found that in the case of SiC>2 the inclusion 
of the C3-terms leads to a significant change of the con- 
tribution involving the off-diagonal terms. In the c^-case 
these contributions are negative, thus they drive the tran- 
sition to lower temperatures. If the C3-terms are taken 
into account, these contributions become less negative 
and hence the transition takes place already at higher 
temperatures. 

In summary, the present calculation shows that MCT 
is able to give an accurate quantitative description of the 



NEP for a strong glass- former. This opens the way for a 
detailed description of the dynamics in the MCT region 
for this class of materials. Hence our results support 
the possibility that MCT is able to make quantitative 
predictions for all types of glass-formers. 
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